/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
/* Copyright (C) 2011 LinBox
 * Written by Brice Boyer <brice.boyer@imag.fr>
 *
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#ifndef __LINBOX_integer_tools_INL
#define __LINBOX_integer_tools_INL

#include "linbox/field/modular.h"
#include "linbox/matrix/blas-matrix.h"
#include "linbox/algorithms/blas-domain.h"
#include "linbox/integer.h"
#include "linbox/field/PID-integer.h"

#include "linbox/algorithms/cra-domain.h"
#include "linbox/algorithms/cra-full-multip-fixed.h"

#include "linbox/randiter/random-integer.h"
#include "linbox/randiter/random-prime.h"

/*******************************/
/*   Matrix/Vector utilities   */
/*******************************/

namespace LinBox
{
	/*! Get the size (in bits) of the biggest element in an integer matrix.
	 * @param sz the size
	 * @param M the matrix
	 */
	template<class IMatrix>
	size_t & MaxElementSize(size_t & sz, const IMatrix & M)
	{
		sz = 0  ;
		size_t tmp_sz = 0 ;
		typedef typename IMatrix::ConstRawIterator Iter ;
		for (Iter Mit = M.rawBegin() ; Mit != M.rawEnd() ; ++Mit){
			tmp_sz = (*Mit).bitsize();
			if (tmp_sz>sz) sz = tmp_sz ;
		}
		return sz ;

	}

	/*! Get the total size (in bits) of the elements in an integer matrix.
	 * @param sz the size
	 * @param M the matrix
	 */
	template<class IMatrix>
	size_t & SumElementSize(size_t & sz, const IMatrix & M)
	{
		sz = 0  ;
		typedef typename IMatrix::ConstRawIterator Iter ;
		for (Iter Mit = M.rawBegin() ; Mit != M.rawEnd() ; ++Mit){
			sz += (*Mit).bitsize();
		}
		return sz ;

	}

	/*! Get the absolute value of the biggest element in an Integer matrix.
	 * @param big the result
	 * @param M the matrix
	 */
	template<class IMatrix>
	Integer & MaxElementAbs(Integer & big, const IMatrix & M)
	{
		Integer pti = 0 ;
		big = 0 ;
		typedef typename IMatrix::ConstRawIterator Iter ;
		for (Iter Mit = M.rawBegin() ; Mit != M.rawEnd() ; ++Mit){
			(*Mit>big)?(big = *Mit):(NULL) ;
			(*Mit<pti)?(pti = *Mit):(NULL) ;
		}
		return std::max(-pti,big) ;

	}

	/*! Get the value of the biggest element in an Integer matrix.
	 * @param big the result
	 * @param M the matrix
	 */
	template<class IMatrix>
	Integer & MaxElement(Integer & big, const IMatrix & M)
	{
		typedef typename IMatrix::ConstRawIterator Iter ;
		Iter Mit = M.rawBegin() ;
		if (Mit == M.rawEnd()) return big = 0 ;
		big = *Mit ;
		++Mit ;
		for ( ; Mit != M.rawEnd() ; ++Mit){
			(*Mit>big)?(big = *Mit):(NULL) ;
		}
		return big ;

	}

}

/**************************************/
/*   Matrix - Vector multiplication   */
/**************************************/

namespace LinBox
{

}

/**************************************/
/*   Matrix - Matrix multiplication   */
/**************************************/

namespace LinBox
{
	// struct FAM_TYPE {
	// _axpy = 10001,
	// _axmy = 10002,
	// _maxpy= 10003
	// };
	// struct FAM_TYPE {} ;
	// struct _axpy : public FAM_TYPE {} ;
	// struct _axmy : public FAM_TYPE {} ;
	// struct _maxpy : public FAM_TYPE {} ;
	// static int _axpy  =1 ;
	// static int _axmy  =2 ;
	// static int _maxpy =4 ;
#define _axmy  1
#define _axpy  2
#define _maxpy 4

	template<int FAM_TYPE>
	struct IntegerCraMatMul {


		typedef Modular<double>     Field;
		typedef Field::Element      Element;
		typedef BlasMatrix<Element> ModularMatrix ;
		typedef BlasMatrix<Integer> IntegerMatrix ;

#ifdef _LB_MM_TIMING
#ifdef _OPENMP
		typedef LinBox::OMPTimer Mytime;
#else
		typedef LinBox::Timer    Mytime;
#endif
#endif

		const IntegerMatrix &_A_, &_B_;

#ifdef _LB_MM_TIMING
		mutable Mytime chrono;
#endif

		IntegerCraMatMul(const IntegerMatrix& A, const IntegerMatrix& B) :
			_A_(A), _B_(B)
		{
#ifdef _LB_MM_TIMING
			chrono.clear();
#endif
			linbox_check(A.getPointer() == _A_.getPointer());
		}

		IntegerCraMatMul(IntegerMatrix& A, IntegerMatrix& B) :
			_A_(A), _B_(B)
		{
#ifdef _LB_MM_TIMING
			chrono.clear();
#endif
			linbox_check(A.getPointer() == _A_.getPointer());
		}

		ModularMatrix& operator()(ModularMatrix& Cp, const Field& F) const
		{
			BlasMatrixDomain<Field>   BMD(F);

			/*  intialisation */
			ModularMatrix Cpp(_A_.rowdim(),_B_.coldim());
			Cp = Cpp ;
			ModularMatrix Ap(_A_, F);
			ModularMatrix Bp(_B_, F);

			/*  multiplication mod p */

#ifdef _LB_MM_TIMING
			Mytime matmul; matmul.clear(); matmul.start();
#endif
			// BMD.mul( static_cast<BlasMatrix<double>&>(Cp),Ap,Bp);
			if (FAM_TYPE == _axpy)
				BMD.axpyin(Cp,Ap,Bp);
			else if (FAM_TYPE == _axmy)
				BMD.axmyin(Cp,Ap,Bp);
			else if (FAM_TYPE == _maxpy)
				BMD.maxpyin(Cp,Ap,Bp);
#ifdef _LB_MM_TIMING
			matmul.stop();
			this->chrono+=matmul;
#endif
#if 0
			if (Ap.rowdim() <= 20 && Ap.coldim() <= 20) {
				Integer chara;
				F.characteristic(chara);
				F.write(cout) << endl;
				cout << "p:=" << chara << ';' << std::endl;
				A.write(cout<< "A:=",true) << ';' << std::endl;
				Ap.write(cout << "Ap:=", F, true) << ';' << endl;
				Bp.write(cout << "Bp:=", F, true) << ';' << endl;
				Cp.write(cout<< "Cp:=", F, true) << ';' << endl;
			}
#endif
			return Cp;
		}

		// BlasMatrix<double> initialize()
		// {
		// return BlasMatrix<double>(A.rowdim(), B.coldim());
		// }


	};

	// template<class Element>
	// struct CRATemporaryVectorTrait<IntegerCraMatMul<enum FAM_TYPE>, Element>
	// {
	// typedef BlasMatrix<double> Type_t;
	// };


	// C = -C
	BlasMatrix<Integer> & NegIn(BlasMatrix<Integer> & C)
	{
		BlasMatrix<Integer>::RawIterator it = C.rawBegin();
		for ( ; it != C.rawEnd() ; ++it)
			Integer::negin(*it);
		return C;
	}

	// D = -C
	BlasMatrix<Integer> & Neg(BlasMatrix<Integer>       & D,
				  const BlasMatrix<Integer> & C)
	{
		BlasMatrix<Integer>::RawIterator      it = D.rawBegin() ;
		BlasMatrix<Integer>::ConstRawIterator jt = C.rawBegin() ;
		for (; it != D.rawEnd() ; ++jt, ++it)
			Integer::neg(*it,*jt);
		return D;
	}

	// C = 0
	BlasMatrix<Integer> & InitZero(BlasMatrix<Integer> & C)
	{
		BlasMatrix<Integer>::RawIterator it = C.rawBegin();
		const Integer zero (0);
		for (; it != C.rawEnd() ; ++it)
			*it = zero;
		return C;
	}

	// C *= b
	template<class T>
	BlasMatrix<Integer> & mulin(BlasMatrix<Integer> & C,
				     const T             & b)
	{
		BlasMatrix<Integer>::RawIterator it = C.rawBegin();
		for ( ; it != C.rawEnd() ; ++it)
			Integer::mulin(*it,b);
		return C;
	}

	// D = b * C
	template<class T>
	BlasMatrix<Integer> & mul(BlasMatrix<Integer>       & D,
				   const T                   & b,
				   const BlasMatrix<Integer> & C)
	{
		BlasMatrix<Integer>::RawIterator      it = D.rawBegin() ;
		BlasMatrix<Integer>::ConstRawIterator jt = C.rawBegin() ;
		for ( ;  it != D.rawEnd() ; ++it, ++jt)
			Integer::mul(*it,*jt,b);
		return D;
	}

	// C += b D
	template<class T>
	BlasMatrix<Integer> & axpyin(BlasMatrix<Integer>       & C,
				     const T                   & b,
				     const BlasMatrix<Integer> & D)
	{
		BlasMatrix<Integer>::RawIterator it =  C.rawBegin() ;
		BlasMatrix<Integer>::ConstRawIterator jt = D.rawBegin() ;
		for (; it != C.rawEnd() ; ++it, ++jt)
			Integer::axpyin(*it,b,*jt);
		return C;
	}

	// D = b E + C
	template<class T>
	BlasMatrix<Integer> & axpy(BlasMatrix<Integer>       & D,
				   const T                   & b,
				   const BlasMatrix<Integer> & E,
				   const BlasMatrix<Integer> & C)
	{
		BlasMatrix<Integer>::RawIterator      it = D.rawBegin();
		BlasMatrix<Integer>::ConstRawIterator jt = E.rawBegin() ;
		BlasMatrix<Integer>::ConstRawIterator kt = C.rawBegin() ;
		for (; it != D.rawEnd() ; ++it, ++jt, ++kt) {
			Integer::axpy(*it,*jt,b,*kt);
			linbox_check(*it == *jt*b+*kt);
		}
		return D;
	}

	// C -= b D
	template<class T>
	BlasMatrix<Integer> & maxpyin(BlasMatrix<Integer>       & C,
					const T                   & b,
					const BlasMatrix<Integer> & D)
	{
		BlasMatrix<Integer>::RawIterator      it = C.rawBegin()  ;
		BlasMatrix<Integer>::ConstRawIterator jt = D.rawBegin() ;
		for ( ; it != C.rawEnd() ; ++it, ++jt)
			Integer::axmyin(*it,b,*jt);
		return C;
	}

	// D = b E - C
	template<class T>
	BlasMatrix<Integer> & axmy(BlasMatrix<Integer>       & D,
				      const T                   & b,
				      const BlasMatrix<Integer> & E,
				      const BlasMatrix<Integer> & C)
	{
		BlasMatrix<Integer>::RawIterator      it = D.rawBegin();
		BlasMatrix<Integer>::ConstRawIterator jt = E.rawBegin() ;
		BlasMatrix<Integer>::ConstRawIterator kt = C.rawBegin() ;
		for (; it != D.rawEnd() ; ++it, ++jt, ++kt)
			Integer::axmy(*it,*jt,b,*kt);
		return D;
	}


	// C += A B
	template<int FAM_TYPE>
	BlasMatrix<Integer> & MMmul_naif_accin_(BlasMatrix<Integer>      & C,
					       const BlasMatrix<Integer> & A,
					       const BlasMatrix<Integer> & B)
	{
		int lda = A.getStride() ;
		int ldb = B.getStride() ;
		int ldc = C.getStride() ;
		// Integer zero(0);
		Integer * A_it = A.getPointer();
		Integer * B_it = B.getPointer();
		Integer * C_it = C.getWritePointer();
		size_t i = 0 ;

		size_t M = A.rowdim();
		size_t K = A.coldim();
		size_t N = C.coldim();

#ifdef _LB_INT_MM_USE_TRANSPOSE
		/*  transpose B */
		IntegerBlasMatrix BT(N,K) ;
		Integer * BT_it = BT.getWritePointer();
#ifdef _OPENMP
#pragma omp parallel for
#endif
		for (i = 0 ; i < N ; ++i)
			for(size_t j = 0 ; j < K ; ++j)
				*(BT_it+i*K+j) =  *(B_it+j*ldb+i);

		/*  mult !!  */
		for (i = 0 ; i < M ; ++i) {
			Integer * C_i = C_it + i*ldc ;
			const Integer * A_i = A_it + i*lda ;
#ifdef _OPENMP
#pragma omp parallel for
#endif
			for (size_t j = 0; j < N; ++j) {
				const Integer * BT_j = BT_it + j*K ;
				for (size_t k = 0 ; k < K ; ++k) {
					if (FAM_TYPE == _axpy)       /*  C = AB + C */
						Integer::axpyin(*(C_i+j), *(A_i+k), *(BT_j+k));
					// else if (FAM_TYPE == _axmy)  /*  C = AB - C */
						// Integer::axmyin(*(C_i+j), *(A_i+k), *(BT_j+k));
					else if (FAM_TYPE == _maxpy) /*  C = -AB + C  */
						Integer::maxpyin(*(C_i+j), *(A_i+k), *(BT_j+k));
				}
			}
		}
#else /*  no pre-transpose */
		/*  mult !!  */
		for (i = 0 ; i < M ; ++i) {
			// Integer * C_i = C_it + i*ldc ;
			// const Integer * A_i = A_it + i*lda ;
#ifdef _OPENMP
#pragma omp parallel for
#endif
			for (size_t j = 0; j < N; ++j) {
				for (size_t k = 0 ; k < K ; ++k) {
					if (FAM_TYPE == _axpy) {       /*  C = AB + C */
						Integer::axpyin(*(C_it+i*ldc+j), *(A_it+i*lda+k), *(B_it+k*ldb+j));
					}
					// else if (FAM_TYPE == _axmy) { /*  C = AB - C */
						// Integer cc = *(C_it+i*ldc+j) ;
						// Integer::axmyin(*(C_it+i*ldc+j), *(A_it+i*lda+k), *(B_it+k*ldb+j));
						// linbox_check(  *(C_it+i*ldc+j) ==  *(A_it+i*lda+k) * *(B_it+k*ldb+j) - cc);

					// }
					else if (FAM_TYPE == _maxpy) { /*  C = -AB + C */
						Integer::maxpyin(*(C_it+i*ldc+j), *(A_it+i*lda+k), *(B_it+k*ldb+j));
					}
				}
			}
		}

#endif

		return C ;
	}

	// C = alpha A * B + beta C
	BlasMatrix<Integer> & MMmul_naif_in( BlasMatrix<Integer>       & C,
					     const Integer             & alpha,
					     const BlasMatrix<Integer> & A,
					     const BlasMatrix<Integer> & B,
					     const Integer             & beta)
	{
		size_t M = C.rowdim() ;
		size_t N = C.coldim() ;
		size_t K = B.rowdim();

		if (!K) {
			if (beta == 0)
				return InitZero(C);
			if (beta == -1)
				return NegIn(C);
			if (beta != 1)
				return mulin(C,beta);
			return C;
		}
		if (!M || !N)
			return C ;

		linbox_check(M == A.rowdim());
		linbox_check(N == B.coldim());
		linbox_check(K == A.coldim());

		/*  non trivial matrices */
		/*  alpha = 0 */
		if (alpha == 0) {
			if (beta == 0)
				return InitZero(C);
			if (beta == 1)
				return C ;
			if (beta == -1)
				return NegIn(C);
			return mulin(C,beta);
		}
		/*  alpha != 0 */
		/*  beta = 0 */
		if (beta == 0) {
			if (alpha == 1) {/*  C = AB */
				InitZero(C);
				MMmul_naif_accin_<_axpy>(C,A,B) ;
				return C;
			}
			if (alpha == -1){ /*  C = -AB */
				InitZero(C);
				return MMmul_naif_accin_<_maxpy>(C,A,B) ;
			}
			else {/*  C = alpha AB */
				InitZero(C);
				MMmul_naif_accin_<_axpy>(C,A,B) ;
				return mulin(C,alpha);
			}
		}
		/*  beta = 1 */
		if (beta == 1) {
			if (alpha == 1) { /*  C =  AB + C */
				return MMmul_naif_accin_<_axpy>(C,A,B);
			}
			if (alpha == -1) { /*  C = - AB + C */
				return MMmul_naif_accin_<_maxpy>(C,A,B);
			}
			/*  C = alpha AB + C */
			BlasMatrix<Integer> D(M,N);
			MMmul_naif_accin_<_axpy>(D,A,B);
			return axpyin(C,alpha,D);
		}
		/*  beta = -1 */
		if (beta == -1) {
			if (alpha == 1) {/*  C =  AB - C */
				MMmul_naif_accin_<_maxpy>(C,A,B);
				return NegIn(C);

			}
			if (alpha == -1) {/*  C = -AB -C */
				MMmul_naif_accin_<_axpy>(C,A,B);
				return NegIn(C);
			}
			/*  C = alpha AB - C */
			BlasMatrix<Integer> D(M,N);
			MMmul_naif_accin_<_axpy>(D,A,B);
			return maxpyin(C,alpha,D);
		}
		/*  beta general */
		/*  alpha = 1 */
		if (alpha == 1) {/*  C =  AB + beta C */
			mulin(C,beta);
			return MMmul_naif_accin_<_axpy>(C,A,B);
		}
		/*  alpha = -1 */
		if (alpha == -1) {/*  C = -AB + beta C */
			mulin(C,beta);
			return MMmul_naif_accin_<_maxpy>(C,A,B);
		}
		/*  alpha general */
		/*  C = alpha AB + beta C */
		{
			mulin(C,beta);
			BlasMatrix<Integer> D(M,N);
			MMmul_naif_accin_<_axpy>(D,A,B);
			return axpyin(C,alpha,D);
		}
	}

	// D = alpha A * B + beta C
	BlasMatrix<Integer> & MMmul_naif( BlasMatrix<Integer>         & D
					  , const Integer             & alpha
					  , const BlasMatrix<Integer> & A
					  , const BlasMatrix<Integer> & B
					  , const Integer             & beta
					  , const BlasMatrix<Integer> & C)
	{
		size_t M = C.rowdim() ;
		size_t N = C.coldim() ;
		size_t K = B.rowdim();

		linbox_check(M == A.rowdim());
		linbox_check(N == B.coldim());
		linbox_check(K == A.coldim());
		linbox_check(M == D.rowdim());
		linbox_check(N == D.coldim());

		if (!K) { /*  A, B are empty */
			if (beta == 0)
				return InitZero(D);
			if (beta == -1)
				return Neg(D,C);
			if (beta != 1)
				return mul(D,beta,C);
			return D;
		}
		if (!M || !N)
			return D ;

		/*  non trivial matrices */
		/*  alpha = 0 */
		if (alpha == 0) {
			if (beta == 0)
				return InitZero(D);
			if (beta == 1)
				return D=C ;
			if (beta == -1)
				return Neg(D,C);
			return mul(D,beta,C);
		}
		/*  alpha != 0 */
		/*  beta = 0 */
		if (beta == 0) {
			if (alpha == 1) {/*  D = AB */
				InitZero(D);
				return MMmul_naif_accin_<_axpy>(D,A,B) ;
			}
			if (alpha == -1){ /*  C = -AB */
				InitZero(D);
				return MMmul_naif_accin_<_maxpy>(D,A,B) ;
			}
			else {/*  C = alpha AB */
				InitZero(D);
				MMmul_naif_accin_<_axpy>(D,A,B) ;
				return mulin(D,alpha);
			}
		}
		/*  beta = 1 */
		if (beta == 1) {
			if (alpha == 1) { /*  D =  AB + C */
				D = C ; //! @todo ou alors \c MMmul_naif_in_<_axpy>(D,A,B,C); ??
				return MMmul_naif_accin_<_axpy>(D,A,B);
			}
			if (alpha == -1) {/*  D = - AB + C */
				D = C ;
				return MMmul_naif_accin_<_maxpy>(D,A,B);
			}
			/*  D = alpha AB + C */
			BlasMatrix<Integer> E(M,N);
			MMmul_naif_accin_<_axpy>(E,A,B);
			return axpy(D,alpha,E,C);
		}
		/*  beta = -1 */
		if (beta == -1) {
			if (alpha == 1){ /*  D =  AB - C */
				Neg(D,C);
				return MMmul_naif_accin_<_axpy>(D,A,B);
			}
			if (alpha == -1) {/*  D = -AB -C */
				Neg(D,C);
				return MMmul_naif_accin_<_maxpy>(D,A,B);
			}
			/*  D = alpha AB - C */
			BlasMatrix<Integer> E(M,N);
			MMmul_naif_accin_<_axpy>(E,A,B);
			return axmy(D,alpha,E,C);
		}
		/*  beta general */
		/*  alpha = 1 */
		if (alpha == 1) {/*  D =  AB + beta C */
			mul(D,beta,C);
			return MMmul_naif_accin_<_axpy>(D,A,B);
		}
		/*  alpha = -1 */
		if (alpha == -1) {/*  C = -AB + beta C */
			mul(D,beta,C);
			return MMmul_naif_accin_<_maxpy>(D,A,B);
		}
		/*  alpha general */
		/*  D = alpha AB + beta C */
		{
			mul(D,beta,C);
			BlasMatrix<Integer> E(M,N);
			MMmul_naif_accin_<_axpy>(E,A,B);
			return axpyin(D,alpha,E);
		}
	}

	// C += A B
	template<int FAM_TYPE>
	BlasMatrix<Integer> & MMmul_cra_accin_(BlasMatrix<Integer>      & C,
					      const BlasMatrix<Integer> & A,
					      const BlasMatrix<Integer> & B,
					      const double              logC)
	{

		// double UpperBound = (naturallog(mA) + naturallog(mB) + naturallog(A.coldim()))+1;

		size_t PrimeSize = 22; //! @todo pourqoi ?

		// IgnobleM = A.rowdim();
		// IgnobleN = B.coldim();

		typedef Modular<double> ModularField ;

		// IntegerBlasMatrix C(A.rowdim(),B.coldim());
		{
			// typedef Modd					RNSField;
			// typedef BlasMatrix<RNSField::Element>		BlasMat;

			RandomPrimeIterator genprime( PrimeSize );
			ChineseRemainder< FullMultipBlasMatCRA< ModularField > > cra( std::pair<size_t,double>(C.rowdim()*C.coldim(), logC) );
			IntegerCraMatMul<FAM_TYPE> iteration(A,B);

			cra(C, iteration, genprime);

#ifdef _LB_DEBUG
			std::cout << "Sole modular matrix multiplications: " << iteration.chrono << std::endl;

			Integer mC; MaxElement(mC, C);
			std::cout << "C max: " << logtwo(mC) <<  " (" << naturallog(mC) << ')' << std::endl;
#endif

		}

		return C;
	}

	// C = alpha A * B + beta C
	BlasMatrix<Integer> & MMmul_cra_in( BlasMatrix<Integer>       & C,
					    const Integer             & alpha,
					    const BlasMatrix<Integer> & A,
					    const BlasMatrix<Integer> & B,
					    const Integer             & beta);

}

/************************/
/*  P-Adic conversions  */
/************************/

namespace LinBox
{

}

#endif // __LINBOX_integer_tools_INL

